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ABSTRACT 

We have developed a new formalism to compute the thermodynamic coef- 
ficient I\ in the theory of stellar and atmospheric stability. We generalize the 
classical derivation of the first adiabatic index, which is based on the assumption 
of thermal ionization and equilibrium between gas and radiation temperature, 
towards an expression which incorporates photo-ionization due to radiation with 
a temperature T ra d different from the local kinetic gas temperature. Our formal- 
ism considers the important non-LTE conditions in the extended atmospheres 
of supergiant stars. An application to the Kurucz grid of cool supergiant atmo- 
spheres demonstrates that models with T ra d— T e s between 6500 K and 7500 K 
become most unstable against dynamic perturbations, according to Ledoux' sta- 
bility integral <Ti>. This results from Ti and <I\> acquiring very low values, 
below 4/3, throughout the entire stellar atmosphere, which causes very high 
gas compression ratios around these effective temperatures. Based on detailed 
NLTE-calculations, we discuss atmospheric instability of pulsating massive yel- 
low supergiants, like the hypergiant p Cas (Ia + ), which exist in the extension of 
the Cepheid instability strip, near the Eddington luminosity limit. 

Subject headings: instabilities — stars: atmospheres — stellar dynamics — pul- 
sations — supergiants — stars: variables: hypergiants — Cepheids 



1. Introduction 

The atmospheres of cool massive supergiants are unstable, which causes pulsation- 
variability, strongly developed large-scale atmospheric motion fields, excessive mass-loss, and 
extended circumstellar envelopes. One of the best studied examples of these very luminous 



- 2 - 



supergiants is the yellow hypergiant p Cas (F2— G Ia + ). This evolved star exhibits stable 
pulsation (quasi-) periods of 300—500 d. 

Although the k- and 7-mechanisms have been identified as the main cause for driving 
pulsations of the less luminous high-gravity atmospheres of Cepheids, little is known about 
the efficiency of these effects for the much more extended and tenuous atmospheres of cool 
massive supergiants. In detailed calculations of the first generalized adiabatic index F\, 
Lobel et al. (1992) (paper I) found that this quantity assumes very small values, below 4/3, 
in low-gravity model atmospheres with 5000 K < T c g < 8000 K, primarily due to the partial 
thermal ionization of hydrogen. Stothers & Chin (1999) recently suggested that enhanced 
mass-loss due to ionization-induced dynamical instability of the outer envelope of luminous 
supergiants which evolve redwards, would terminate their redward movement, and provide 
an explanation for an observational lack of yellow and red supergiants with log(L*/L Q )>6.0. 

However, the major problem for evaluating supergiant dynamic (in)stability, based on 
detailed calculations of Ti, is the breakdown of LTE conditions in these very extended 
atmospheres. The importance of NLTE ionization- and excitation-conditions is evident from 
modeling the spectra of these stars, which are formed in conditions of very small gravity 
acceleration. The local ionization equilibrium is strongly determined by the stellar radiation 
field, which determines important thermodynamic quantities such as the heat capacities, and 
the related mechanic compressibility of these atmospheres. 

In this paper we develop, for the first time, a self-consistent thermodynamic formalism 
which accounts for departures from LTE for the calculation of lY This goal is accomplished 
by introducing the temperature of the radiation field as an independent state variable which 
can differ from the local kinetic gas temperature. We discuss in Sect. 2 the departure 
from thermal ionization equilibrium by an incident and diluted stellar radiation field in the 
Eddington approximation. Section 3 provides a historical overview of the development of the 
theory of the adiabatic indices. The complete analytical expressions for the computation of 
Ti and the heat capacities in mixtures of monatomic gas, interacting with radiation are given 
in Sect. 4. We account for departures from LTE conditions due to the interaction of matter 
and radiation, by evaluating the thermodynamic quantities accordingly. Section 5 presents 
a discussion of the effects of NLTE conditions on r x in cool supergiants. These detailed 
NLTE calculations of I\ are applied in Sect. 7, to evaluate their dynamic stability according 
Ledoux' stability integral <Ti> for radial fundamental mode oscillations (Sect. 6). We 
apply our calculations to a new (Kurucz) grid of cool supergiant model atmospheres, which 
we compute down into the stellar envelope. We will demonstrate that NLTE-ionization of 
hydrogen strongly enhances the destabilization of supergiant atmospheres with 6500 K < T c s 
< 7500 K. For models towards smaller gravity acceleration the stability integral decreases, 
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and the destabilizing regions occur at lower densities over a larger geometric fraction of 
the atmosphere. A discussion of these results in relation to pulsation driving in yellow 
hypergiants, and atmospheric instability regions recently identified in the upper portion of 
the HR-diagram by de Jager & Nieuwenhuijzen (1997), is given in Sect. 8. The conclusions 
of this theory and application are listed in Sect. 9. 



In the atmospheres of supergiants the ionization balance of tenuous gas is not solely 
determined by a collisional equilibrium according the Saha-equation. Significant departures 
from thermal (equilibrium) ionization occur due to photo-ionization by an incident radiation 
field. The local ionization state then becomes dependent of both kinetic gas temperature, 
and the temperature of the radiation field. For example, stellar UV radiation strongly 
influences the Balmer continuum in hot stars. In deeper atmospheric layers, where the 
atmosphere is sufficiently optically thick for all wavelengths, both temperatures thermalize, 
and the equilibrium radiation field assumes an intensity distribution determined by the local 
kinetic gas temperature. In the upper layers the radiation field dilutes with distance from a 
point in the atmosphere where the gas becomes sufficiently optically thin, and the radiation 
temperature decouples from the local thermodynamic conditions. 



A comprehensive description of partially ionized systems which deviate from equilibrium, 
due to a radiation field of temperature T ra d, not in equilibrium with the electron Maxwell 
distribution of temperature T e , is given in Elwert (1952). All particle components (neutrals, 
ions and electrons) are assumed to be in a Maxwell distribution. For the calculation of ion- 
ization fractions, this statistical theory assumes that the detailed balance between collisional 
and/or radiative ionization, and recombination processes applies. It enables to express the 
ionization fractions through a departure coefficient b from the Saha equilibrium (also the 
'NLTE Saha-equation'), which can be evaluated using a reduced form of the collision ion- 
ization cross-section and, for the photo-ionization coefficient, a diluted Planck distribution 
with a cross-section obtained from quantum-mechanical calculations. The Elwert-equation 
is given by: 



2. 



Non-LTE ionization 



2.1. Eddington approximation 




(1) 
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where rij denotes the number density of particles in the j-th ionization stage, and n e is the 
electron number density. The Saha-Boltzmann equation for thermal ionization form the r-th 
excitation level is: 
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n j+in B \ _ 2 uj+i ffrj+i / 27rm e kT e y ^ / Ij - Xr, 3 , ^ 
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with the partition function Wj = S^j ex P ( — Fr") > anc ^ h ^ ne ionization energy from the 

ground state. Xr,j is the excitation energy of level r and g r j its statistical weight, h is the 
Planck constant and the other symbols have their usual meaning. The departure coefficient 
(for ionization from the ground level) in Eq. (1) is: 



bj = 1 * , (3) 



1 + W -£r c t ex P(l/c - Z/rad) ^(j/e, Z/rad) 

with the function: 

fl--) 

F(y c ,y rad )^ \ ^ , (4) 

I 1 "ij 

where we denote; y e = Ij/ (kT e ) and y ra( j = ij/ (A;T ra d). A and £> are constants which depend 
on the ionization energy, the Thompson cross-section, and the Bohr radius. Note that Eq. 
(2) depends only on the kinetic gas temperature T e , whereas Eq. (1) is dependent of T ra( j as 
well. The dilution factor with geometric height d from the stellar surface R+ is: 



W(z) = 1/2 (1 - y/1 - 1/z*) , (5) 

where z=d/R ic . At the surface z—1, the dilution factor W(z=l) = l/2, because a gas particle 
is irradiated at most by half the stellar hemisphere. 

Ecker (1978) distinguished two important conditions of partial ionization from the gen- 
eral equation (1), based on a critical electron density n c , which is a function of the kinetic 
temperature, and the ratio of the radiation and kinetic temperature: 

i 

B yl 

n c (T c , T rad /T c ) = — — exp(y c - y rad ) F(y e , y rad ) . (6) 

A 2/rad 

%. The Corona case assumes that the radiation density is so small that photo-ionization 
is negligible in comparison with electron collision-ionization, and the electron density is still 
so small that three-body recombination is negligible in comparison with radiative recombi- 
nation: 

n c (T c , 1) > n c > Wn c (T c ,T rad /T e ) . (7) 
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Hence the departure coefficient can be approximated by: 

which demonstrates that for coronal conditions the plasma becomes 'under '-ionized to a 
degree smaller than the Saha-equilibrium, and for which the ionization fraction becomes 
independent of the local electron number density, since it cancels out in Eq. (1) with Eq. 
(8). 

%%. The Eddington case assumes that the electron density is so small that ionization 
is dominated by photo-ionization, and three-body recombination is dominated by radiative 
recombination: 

n c (T c ,l)>n c and W n c (T c ,T rad /T c ) > n e . (9) 
Hence, the departure coefficient can be approximated by: 



1 1/rad ( , F(y e ,y e ) 

— exp(y rad - y e ) — 

W y c F(y e ,y Tad ) 



h ^ ^ y —ex P (y rad - y e ) -^J±L (i ) 



! Ijl (k (± l _W v 1 "^/ i • 

" WT TaA 6XP ^ k W d tJ) (i _ ■ 1 j 

Note that for these conditions of relatively low electron density and high radiation den- 
sity, the ionization fraction still depends on the electron temperature through the radiative 
recombination mechanism. 



2.2. Governing ionization equation 

For our calculation of T 1 in the atmospheres of supergiants we consider the Eddington 
approximation, for which photo-ionization dominates collisional ionization. The departures 
from LTE become very large in low density, optically thin, regions. Because of the radial 
variation of the local temperature and the wavelength variation of the opacity sources, the 
emitted spectrum departs from a blackbody at any particular temperature. However, for 
T e ff above 4000 K, the local gas density and temperature are only weakly determined by 
the ambient radiation field, because important molecular opacity sources remain limited 
for these conditions. For our calculations we can assume that the stellar radiation field 
in the atmosphere where tr oss <2/3 (i.e. above R+), can be approximated by the effective 
temperature; T ra( j— T e g. For 'grey' atmospheres in radiative equilibrium, the color (surface 
brightness or radiation) temperature is 0.811 x T cS (Woolley & Stibbs 1953). 
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We also consider non-LTE ionization conditions in atmospheric regions where tr oss <2/3, 
with T c and T rad below 20,000 K, for ionization energies Ij in excess of 7 eV. Hence, the 
trailing factor at the right-hand side of Eq. (11) approaches unity or: 

6 ^±A exp (^M IV). (12) 

It represents a more tractable expression for our further derivation of important thermody- 
namic derivatives in Sect. 4.2. The NLTE-Saha equation in the Eddington approximation 
is hence obtained from Eq. (1), Eq. (2), and Eq. (12), for single ionization (j—0— 1) from 
the ground level (or Xr,j — 0) of element i: 



n ~ u ,i \ h 2 J \ kT iad J \T, 



'rad 



Eddington (1926) first considered ionization of the interstellar medium due to ionizing 
radiation from nearby stars, and obtained: 
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n U ,i \ h 2 J V kT rad 

with T rad of the order of the effective temperatures of stellar atmospheres. The factor at the 

i 

right-hand side of Eq. (13), (T e /T rad ) 2 is due to thermal motions, and first appears in Rosse- 
land (1936). Stromgren (1948) gave a further refinement of Eq. (13) for photo-ionization 
from the ground state, which includes recombinations onto energy levels above the ground 
level. Weymann (1962) applied Eq. (13) in a study of the ionization equilibrium in the upper 
atmosphere of the supergiant a Ori (M2 lab), where T e <5500 K. Similar applications to the 
wind conditions of stellar chromospheres are for example given in Hartmann & McGregor 
(1980). 

We conclude this section by emphasizing that the condition of detailed balance, for 
photo-ionizations by an equal number of recombinations (per unit volume and unit time) in 
the Eddington approximation, is required for applying Eq. (13). Next to the condition of 
detailed balance, our calculations of thermodynamic quantities also require that the kinetic 
temperatures of the neutrals, ions, and electrons thermalize on a time-scale shorter than 
the characteristic time-scale of heat exchange within the fluid; i.e. due to an atmospheric 
temperature gradient. The proton-electron relaxation-time (Spitzer 1972) is: 



^^ jaww ^^Z? ., (i5) 
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where Z p = 1 the proton charge number, e the electron charge and A p the proton mass m p in 
atomic units (~ 1). For 2xl0 3 K <T C < 2xl0 4 K, in the atmospheres of supergiants, where 
n c ~10 12 cm -3 , and which are mainly composed of hydrogen, we compute for In A ~ 10 (as 
tabulated in Spitzer 1962), that large hydrodynamic perturbations of the local state variables 
should not occur on time scales shorter than 9xl0~ 7 s<t s <5xl0~ 6 s, or on characteristic 
length scales smaller than Va ow x t s . Equilibrium thermodynamic conditions can for ex- 
ample not establish within the thin layer trailing strong shock waves, where the electron 
temperature departs from the heavy particle gas temperature, or by the presence of strong 
electro-magnetic fields which can separate the neutral and charged particle temperatures of 
a tenuous plasma. 



3. Stellar stability coefficients 

Eddington (1918, 1919, 1926) first derived an analytical expression for the first adiabatic 
index F 1 , for a mixture of material gas with radiation. His equation (129.52) is: 

r „_,. (4-3OTT-1) 

T ' =0+ fi+iXi-m-m ' (18) 

where 7 is the specific heat ratio, and (3 the ratio of the material gas pressure to the total 
pressure (using modern symbols). He distinguished a general, or 'effective', ratio of specific 
heats from the exponent to the density in the polytropic equation of state: P = Kp r , where 
K is a constant. Strictly thermodynamically speaking, the quantities V and I\ are not 
identical, but the terminology of 'an exponent' to indicate Eddington's adiabatic quantity: 



r, S (17) 

ad 



dlnp 



has been adopted since in the astrophysical literature. Eddington's consideration of evaluat- 
ing adiabatic thermodynamic derivatives for studying adiabatic stellar oscillations was later 
more rigorously addressed by Chandrasekhar (1939). He defined the two other adiabatic 
indices: 



and 



<91nP 
dh\T 



(18) 



dlnp 



r 3 " 1 = > (19) 



ad 



and also first obtained detailed expressions for a mixture of material gas and radiation. 
P denotes the total pressure, i.e. the sum of partial material gas pressures and radiation 
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pressure. Only two of the three adiabatic indices are independent because a non-degenerate 
gas state is determined by at least two independent state variables: 

r 3 -i^ (r V )ri . (20) 

1 2 

With these definitions, Chandrasekhar also obtained more general expressions for the specific 
heats C p and C v , and clearly distinguished I\ from C p /C v . His equation (148); 

ri = (21) 



where 



and 



Cp = ^ {(3 2 + ( 7 - 1)(4 - 3/3) 2 + 12( 7 - 1)0(1 - /?)) , (22) 

C„ = ! G9 + 12(7-l)(l-/3)) , (23) 

with c p and c„ the specific heats of the material gas, marks an important step towards a 
self-consistent description of thermodynamic derivatives required in the theory of dynamical, 
convective, and pulsational stability of stellar atmospheres. Fowler & Guggenheim (1925) 
were the first to derive expressions for the specific heats with radiation, where various stages 
of ionization of the material are allowed. However, they made certain assumptions about the 
weight factors and the excitation of atoms and ions, which were too restricted. Independently, 
Moglich, Riewe, & Rompe (1939) derived expressions for the specific heats of a singly-ionizing 
one-component monatomic gas without radiation, which were corrected shortly after by 
Biermann (1942). He calculated: 

c P = ^(l+x)(^ + x(l-x)(^ + 1 Lj |, (24) 

and 

* = —(! + *) (3 + __(- + _ ] |, (25) 

where I is the ionization energy, x the ionization fraction, and N is the number of atoms 
per unit mass. Biermann also derived the correct expression for the first adiabatic index 
(evidently, without using this terminology), for negligible radiation pressure: 

r,= 5 + *(i-x) + (26) 



3 + *(!-*) (i + d + jL) 2 ) 
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Note that the expression simplifies to Eq. (16) without partial ionization (x—0), since f3— >1 
with vanishing radiation pressure, and hence 1^— >7 = | for monatomic gas. 

Rosa & Unsold (1948) extended the expressions for the specific heats by considering a 
mixture of partially ionizing hydrogen and helium gas. They also provided a detailed numer- 
ical evaluation of their expressions. These calculations followed an investigation by Unsold 
(1938) of Schwarzschild's convection criterion. He improved upon earlier work by Sieden- 
topf (1933a, 1933b, 1935), and obtained the correct equation for the adiabatic temperature 
derivative: 

r 2 = 5 + x(i-x)(| + ^) 2 
r 2 -i 2 + x(i-x)(| + ^) ' 1 ; 

Unsold (1938) introduced the 'mean' degree of ionization x = Y1T ViXi, with the abundance 
of element i, and hence x = 1 for a fully ionized gas. It enabled to obtain an extended ana- 
lytical expression for in terms of summations over the ionization fractions of m elements 
of a gas mixture [see his equation (93,20)]. In the second edition of his monograph on stellar 
atmospheres, Unsold (1968) also derived an expression for c p , by means of the definition of 
x. A similar treatment for c v is given in Menzel, Bhatnagar, & Sen (1963). These equations 
however omit the important influence of a radiation field, as Chandrasekhar demonstrated. 
This problem has been investigated by Krishna-Swamy (1961), who computed the adiabatic 
temperature derivative for a multi-component mixture of singly-ionizing monatomic gas and 
radiation, using an equation of state of the form: P t = N k T (1 + x) p + | a T 4 . The complete 
analytical expressions for the specific heats of this mixture were first derived by Mihalas 
(1965): 



Nk 



= (l + 20a + 16a 2 ) {l+x) + J2 "M 1 ~ x i) (J^j 



x— < x 2 > 

x- < x 2 > I _ _, /5 \ . L t 

x— — \ x[ ' 



x 2 — 2x— < x 2 > 



(l + x)Q + 4a)-£i/ i s i (l-xO^J > (28) 



and 



Nk 



+ ( X ~ X> ' X ' ] kf 2x- <x 2 > ' (29) 



-10- 



where < x 2 >= Yl% u i x h h ls the ionization energy of element i, a is the ratio of the 
radiation pressure and the material gas pressure, and the other symbols have their usual 
meaning. These equations simplify to Eq. (24) and Eq. (25) in case the radiation pressure 
vanishes (a — > 0), for a one-component gas m—1 (hence x = x and < x 2 >= x 2 ). They 
also simplify to Eq. (22) and Eq. (23) when partial ionization of all elements vanishes 
(xi — > 0), since a = 1/(5 — 1, and for monatomic ideal gas c p = 5/2 N k, c v = 3/2 N k, or 
7 = c p /c v = 5/3. The analytical representation of Eq. (28) and Eq. (29) appears rather 
complex, but we will show in Sect. 4.2 that their further generalization, in which the gas 
temperature differs from the radiation temperature, enables to define for every element i two 
functions Gj and Hi, which reduce cp t and c v to basically two terms. 

The calculation of the adiabatic indices for real gas in the context of equilibrium ther- 
modynamics was first given by Cox & Giuli (1968). The specific heat ratio is related to the 
ratio of the isothermal and adiabatic compressibility coefficient via the Maxwell relations: 

CjL = — , (30) 

with the coefficients: 

where S denotes the thermodynamic entropy function. Hence, I\ can be expressed by: 

d In P t \ ( d In P t 



^ X \ din p ) s c v \ d\np ) T ' 

The thermodynamic derivative ^ jftfp ) ' ^ s ^ ne P~^ c ^ OT a ^ the right- hand- side of Eq. (21), 
which Chandrasekhar distinguished from the specific heat ratio by calculating Eddington's 
index for a mixture of ideal material gas and radiation. Cox & Giuli (1968) proposed to term 
this factor ^ d Q\^p ) the "density exponent in the pressure equation of state" xpi probably 
following the adopted designation of 'adiabatic exponent' for Ti, 



d\nP t \ 



^UwJr' hencer i = ^>- ( 33 ) 

Although this terminology has widely been adopted in the literature, we note that this 
quantity is by no means a true 'exponent' in the polytropic equation of state. It is merely a 
coefficient to the specific heat ratio, required to determine an important adiabatic derivative. 
Fowler & Guggenheim (1925) also calculated this quantity and aptly called it 'the isothermal 
factor' (see their Eq. 9), since it is inversely proportional to the isothermal compressibility 
coefficient: k t = (P t x p ) _1 . Cox & Giuli (1968) also introduced the "temperature exponent" 
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(more adequately 'the isochoric factor'): 

which provides an important general relation for the thermodynamics of real gas: 

Pt XT /ocr\ 

I Xp 

Thermodynamic stability demands that both heat capacities and both compressibilities are 
positive. For every real gas cp t > c v and kt > i^s > 0, and cp t —c v does not equal the universal 
gas constant. The latter equality applies only to simple ideal gas, for which xt — Xp — 1- 

With the definitions of the isothermal and isochoric factor, the two other adiabatic 
indices are obtained with: 

r 2 c Pt 



1 2 - 1 c Pt 

and 



XT, (36) 



r 3 - 1 = C 2^L (37) 

Cv XT 

which yields the general identity Eq. (20). The complete expressions for Xp and Xt for a 
singly-ionizing multi-component monatomic gas with radiation are (Lobel et al. 1992): 

Xp (1 + x) (x + Ei 1/^(1 -x*))' 1 } 

and 

PxY], -ViXAl - xM\ + Sr) , , 

Xt= 4-3/3 + n ,-w-lV n VT - 39 
(1 + x) (x + X^W -»0) 

Note that for a hypothetically non- ionizing gas (xj, x) — >0, whence Xp~ > P an d Xr^4— 3/3. 
These expressions simplify to those of Cox & Giuli (1968) for a gas composed of one ionizing 
element (x = x) and radiation. Note that Biermann (1942) first derived Xp without radiation 
((3 — 1), however expressed as a factor to cpjc v : 

Xp = {~8h^) T = 1 ~ (dh^) T = (l + x)(2-x)' (40) 

where \x is the mean molecular weight. In the partial ionization zone of an abundant element 
the mean molecular weight reduces because \i = /io/(l + x), where /io is the mean molecular 
weight of the unionized gas (e.g. /xo=l-26 for the abundance of cosmic material). In these 
regions the mean ionization fraction x increases (0 < x < 1) and, as can be seen from Eq. 
(38), Xp assumes values below unity. The increased compressibility in the thermal ionization 
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regions reduces with Eq. (32) the value of the first adiabatic index r\ to below the monatomic 
gas value of 5/3. Here, compression energy is mainly converted into ionization energy, which 
also changes the local heat capacities. An equilibrium radiation field reduces (3 — > 0, and 
hence Xp ~ > and xt —> 4. 

Lobel et al. (1992) demonstrated that I\ can assume values below unity when isotropic 
radiation pressure is important in the partial ionization region of an abundant element, al- 
though / ~f=cp t /c v > 1 for every stable gas. The higher compressibility of radiation, compared 
to pure material gas of the same temperature and pressure, diminishes Xp to very low values 
for 6100 K< T e <9000 K, in the partial ionization region of hydrogen. Hence, I\ can decrease 
to very small values of 0.84. Note however that when radiation vanishes, the product of the 
specific heat ratio and x P always exceeds unity (I^M). The extended Eqns. (28), (29), (38), 
and (39) also enable to correctly compute the compressibility of atmospheric regions where 
simultaneous ionizations of various elements (e.g. H — > H + and He — > He + ) considerably 
lower the values of ^1,2,3- 

The combination of ionization and radiation can diminish F 1 to below 4/3, which plays 
an important role in the study of dynamic stability of gas spheres, which we discuss in 
Sect. 6. Conventional calculations of the adiabatic indices assume however that the gas and 
radiation temperature are equal, and that no interaction occurs between the material gas and 
the radiation field. We presently investigate the effects of a radiation field on the adiabatic 
indices due to photo-ionization in the Eddington approximation. The local ionization state is 
no longer determined by pure collisional processes, but by the incident and diluted radiation 
field as well. Radiative deviations from the Saha equilibrium produce important effects on 
the overall compressibility of the plasma, which directly determines the dynamic stability of 
supergiant atmospheres. 



4. First adiabatic index with photo-ionization 

4.1. Definition of the generalized functions 

When the radiation temperature differs from the local kinetic gas temperature, the clas- 
sic equation of state, which assumes thermal equilibrium between material gas and radiation 
P t = N k T (1 + x)p + \ a T 4 , is replaced by: 

W 

P t = NkT c (l + x)p + — aT r \ d , (41) 

where the first term is the material gas pressure P g , and the second term the diluted radiation 
pressure P ra d- The gas state is hence determined by three independent state variables, instead 
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of two. Therefore, the calculation of the first adiabatic index must consider T e and T rad as 
independent state variables. The constituent heat capacities 

Cf ' = (m) p ' mdc " = (w)j (42) 



and 



are replaced by: 

Cft= (J£) + OL) , (43 ) 

where e denotes the internal energy function, and h is the enthalpy function. The lower case 
symbols denote 'specific' quantities, or expressed per unit mass. The isothermal factor 

is replaced by: 

It is important to note that the three adiabatic indices become fully independent by the 
introduction of an additional state variable T rad . Since the second and third adiabatic indices 
are temperature derivatives they are also re-defined by: 

r »« = ( B * P >\ and F - = , (47) 



and by 



r 2g -l \d\nTj ad : r 2rad -l V«91nT rad/ad 

. fd\nT e \ fd\nT iad \ 

With these definitions the identity Eq. (20) no longer applies, because it is only valid for 
thermodynamic systems with a unique temperature. 



4.2. Derivation of the generalized functions 

In the Appendix we obtain the detailed expression for the specific heat capacities for 
which T e differs from T rad . The appendix is self-contained and can be read without fur- 
ther cross-references. All thermodynamic quantities and functions are defined, and their 
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detailed expressions are presented there. Major intermediate results, required to obtain the 
detailed expression for cp t and c v , are also provided. We denote 9 = T c /T rad , and find after 
considerable algebra: 



||= Q + 4« (40(a + l) + l)) {l + x) + Y t u i x i ^-x i )H i , 



and 



c v 

with the functions: 



Q + 12 a (1 + x) + "M 1 ~ x i) °i > 



(49) 



(50) 



X= ^1/^(1 -x^ , (51) 

i 

F= (52) 



and for every element i: 



i/,= Wi + 1 + 4a g) l ^^-L^. + JL.\. (55) 



«'= 5 + + (53) 

(l + Aa)x(l + x) -Y U 
x(l+x)+X + ~kT, 

It can be shown that for 9 = 1, Eq. (49) and Eq. (50) simplify to Eq. (28) and Eq. (29). 
The first term in Eq. (49) and Eq. (50) is the heat capacity due to the translational motion 
of neutral atoms, ions, and electrons, and a is dependent of the diluted radiation pressure. 
The second term increases the heat capacities due to extra internal degrees of freedom, which 
results from partial photo-ionization by the incident radiation field. When T ra( j — > T e , the 
NLTE-ionization balance assumes the thermal Saha-equilibrium, and the generalized heat 
capacities simplify to the heat capacities for a unique temperature T. 

The analytical expression for the isothermal factor Xp, m which T ra( j differs from T e , is 
formally identical with Eq. (38): 

= g (x 2 + X + Y^j V%Xi{l - Xi)) 

*p- (l + x)(x + J2 i u l x l (l-x, l )) ■ [ } 

However, the dilution of the radiation pressure enters this generalized expression through (5. 
The second term at the right-hand-side of Eq. (46) vanishes because P rad is invariable for 
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constant T ra d, and P ra< j is independent of T e . The derivative in the first term is evaluated for 
variable T rad , because the radiation temperature determines the NLTE-ionization fraction. 

In Sect. 7 we investigate numerically the properties of the NLTE T 1 for multi-component 
monatomic gas with radiation. If we consider only one element of material gas (m=l and 
x — x), Eq. (49) simplifies to (see Lobel 1997): 



(1 + x) \(\ ■ \n(\0{n 1)1 



Nk v ' {\2 

+ ^{bw.^{H l+ ^y^>- (57) 



and Eq. (50) to: 



|_.(3 + 12ae ) (1+l) + ^)(3 + _L)g +e ( 1 + _I_)). (58) 
Because the isothermal factor Eq. (56) simplifies to: 

Xp = (l+x)(2-x) ' (59) 
we obtain after some algebra the NLTE 'one-component' first adiabatic index: 

c v 

§ + Aa(A9(a + 1) + 1) + f (1 - x) (f + $ + 4a) (| + 6V + 4afl) 

| + 12afl + f(l - x) (Q + $) (± + M>) + § + 12a0) ' ( ' 

where we denote: 



$ = 1 + T7F and * = 1 + rip- . (61) 

fc J e fti rad 

For T c =T rad =T, hence 0=1 and $ = ^ = 1 + ^, it follows that Eq. (60) simplifies to (e.g. 
Eq. (71.16) of Mihalas & Weibel Mihalas 1984): 

| + 12a + f(l- a ;)((| + ^) 2 + | + 12a) ' 

which simplifies without radiation (a—0 and (3=1) to Eq. (26). 

We have also obtained the detailed expressions for the two other adiabatic indices T 2 
and T 3 along similar lines, defined with Eq. (47) and Eq. (48), but which will be presented 
elsewhere. These generalized functions do not simplify to the classic expressions derived 
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for a unique T, because the partial temperature derivatives, computed with the Saha equi- 
librium equation, differ from those obtained using the NLTE-Saha equation. It should be 
remembered that the latter is an approximate expression for which the condition of detailed 
balance is assumed to sustain the condition of (ionization-) equilibrium. The (LTE) Saha- 
equation, however, is exact and follows from the minimization of the free energy function, 
which dictates the condition of thermodynamic equilibrium. 

Note that for the computation of thermodynamic quantities free energy minimization 
techniques have been applied to include the effects of pressure ionization due to Coulomb 
forces, and the effect of degenerate states (see Hummer & Mihalas 1988; Mihalas, Dappen, 
& Hummer 1988; Dappen et al. 1988). This numerical approach offers the advantage, for 
example, to incorporate the temperature dependency of the (properly truncated) internal 
partition functions, and to compute their effect on the heat capacities and the compress- 
ibilities. However, the free energy minimization scheme is strictly numerical because a large 
system of non-linear equations is to be solved iteratively for astrophysical mixtures (e.g. 
Mihalas et al. 1990). 

In contrast, our formal approach for mixtures with simultaneously ionizing elements 
and T rad ^Te, offers the advantage to track analytically the complex behavior of the adi- 
abatic indices. For example, the remarkable occurrence of conditions where Ti<l. Our 
formalism applies however to supergiant atmospheres, for which electrostatic interactions 
and degenerate states are negligible. Possibly, relativistic or degenerate gas effects on Ti can 
be incorporated by the analytical approximations outlined in Cox & Giuli (1968), Beaudet & 
Tassoul (1971), Elliot & Kosovichev (1998), and Stolzmann & Blocker (2000). Ecker & Kroll 
(1963) presented an interesting statistical method, based on the minimum entropy produc- 
tion principle of irreversible thermodynamics, to compute the lowering of ionization energy 
for a Saha-type equation which considers Coulomb interaction. Note also that Mollikutty, 
Das, & Tandon (1989) computed effects of a uniform magnetic field on the adiabatic indices, 
and provided generalized LTE expressions for the heat capacities, including the magnetic 
pressure. 



We conclude this section with an important thermodynamic relation for IV From the 
first law the relationship among P t , p, and e along an isentrope (dS=0) is: 



4.3. General thermodynamic relation for T 



i 




(63) 
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where h = e + P t / p. The introduction of an additional state variable T rad defines surfaces of 
constant (total) entropy with the specific entropy function s=s(p, T e , T rad ) in the space of 
three independent state variables. For adiabatic changes: de— Pt/ p 2 dp and dh — 1/p dP t , 
the first adiabatic index can be expressed by: 



The value of r\ is related with the slope along an adiabate in the (P t , p)-plane, for a 
given T rad or T e . Since the NLTE-Saha equation considers local interactions of material gas 
with the incident radiation field, the total entropy function cannot classically be obtained 
by summing over the particle entropy fractions, the entropy of released electrons, and the 
entropy of equilibrium radiation (i.e. by means of the Sackur- Tetrode equation; see Lobel 
1997). s(p, T e , T rad ) is to be obtained from the temperature derivative of the free energy, 
defined for the gas state determined by T c and T rad . Similarly, for conditions of equilibrium, 
e=e(p, T e , T ra d) and h—h(p, T c , T rad ) are uniquely defined functions, and it can be shown 
that their ratio 7h = h/e is related to the first adiabatic index by: 



Hence, the decrease of Ti due to partial ionization and the release of bound electrons is 
related with the decrease of 7^ (e.g. Nieuwenhuijzen et al. 1993). While the value of I\ 
is related with the adiabatic compressibility of the fluid (e.g. the bulk modulus) by local 
mechanic perturbations, the value of 7h is determined by the corresponding changes of the 
internal energy. 



Figure 1 shows the behavior of Ti and the mean ionization fraction x, for variable T e and 
^rad- We compute the thermodynamic quantities for a mixture of 16 elements comprising 
H, He, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, K, Ca, Cr, and Fe, for solar abundance values 
(Anders & Grevesse 1989). The internal partition functions are derived according the meth- 
ods developed by Claas (1951) and by Baschek, Holweger, & Traving (1966). The left-hand 
panels display the NLTE Ti computed at the stellar radius (W—l/2), and the right-hand 
panels show x. In the upper panels the material gas pressure is set to 10 dyn cm" 2 , whereas 
0.5 dyn cm" 2 is used for the lower panels. These conditions are typical for the atmospheres of 
cool supergiants. Figure 1 reveals that T x acquires minimum values for T rad between 7000 K 
and 9000 K, primarily due to the partial ionization of hydrogen. Towards very low T rad , Ti 
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5. Behavior of T± in supergiant atmospheres 
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approaches the value of 5/3 for monatomic ideal gas, whereas for high T ra d Ti assumes the 
equilibrium radiation value of 4/3. The intermediate regions with lowest Ti correspond to 
a;~0.5. The NLTE-ionization equation causes a strong dependence of the local mean ioniza- 
tion fraction on the value of T ra( j. For a given radiation temperature, however at constant gas 
pressure, the ionization fraction decreases towards higher kinetic gas temperatures, because 
the NLTE-ionization fraction remains dependent of the local kinetic temperature through 
the factor Ti in Eq. (13) (but which is Ti for pure collisional ionization). 

We find minimum values for I\ which can decrease to below 0.7 for very small gas 
pressures and low T e . The minima gradually increase for larger gas pressures, by shifting 
towards higher T ra( j. This is also shown by the resurface plots of Fig. 2. The deep minimum 
around T rad ~8000 K (lower panel) becomes shallower and broadens for increased kinetic 
pressures (upper panel). The relative decrease of the radiation pressure also enhances the 
partial ionization of He, which is visible for T rad between 14,000 K and 16,000 K. Note that 
we have also included the partial ionization of He + in our calculations. Its ionization can 
be treated as a separate 'element', because the ionization energy is very high, and partial 
ionization occurs around 25,000 K for -P g =l dyn cm~ 2 . Electrons from less abundant metal 
atoms, with smaller ionization energies, contribute only slightly to the total electron pressure. 
The much larger abundance of hydrogen (by about a factor of ten larger than He) causes the 
large decrease of T± to occur for T ra d between 7000 K and 9000 K, for the small gas pressure 
conditions of supergiant atmospheres. 

An important influence on the local mean ionization fraction is the dilution of radiation 
pressure with distance above i?*, which we further discuss in Sect. 7.1. Calculations with 
W(z)<l/2 (for z > 1) reveal that the minima of Ti increase, which is expected when kinetic 
pressure dominates the gas state. 



6. Theory of stellar dynamic stability 

The theory of adiabatic oscillation of gas spheres shows that, for homologous radial 
motion, the square of the angular frequency of the lowest radial pulsation mode, in the 
linear approximation, is given by (Ledoux 1945, 1965): 

a ° = l-R 2 A\r =<ri "3 >- /~' (66) 

J r 2 pdV 6 J o 

where Qo is the gravitational potential of the star in its equilibrium state, and I denotes 
the moment of inertia with respect to the center of the star. Since Ti is a function of depth 
in the star, the integration requires geometrical depths down to the stellar center for the 
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computation of <7q. When <Xq > the configuration is stable, because the standing wave 
solution of the equation of motion does not grow with time. For supergiants we can assume 
that the fundamental eigenfrequency a is not affected by magnetic fields, and that rotational 
kinetic energy can be neglected. 



<Ti>= (67) 



Stellar dynamic stability depends on the volume-averaged value of Ti, and vanishes 
when this average 

is less than 4/3, hence aft < 0. Because our detailed knowledge of Ti is limited to the 
atmospheric layers, it is not straightforward to infer stellar dynamic stability from Eq. (67). 
However, Stothers (1999) recently argued that the integral can be limited to 'an effective 
basis', or a point where <I\> becomes constant in deeper parts of the atmosphere near the 
base of the stellar envelope. Numerical integrations of the equation of motion demonstrate 
that this truncation is allowed, because the radius eigenfunction S(r)/r (the relative pulsation 
amplitude) at the base of the outer envelope is already many orders of magnitude smaller 
than its value at the stellar surface (so small that it is virtually indistinguishable from zero). 
Since regular stellar pulsation is caused by envelope mechanisms, this limitation of Ledoux' 
stability integral provides a useful means to test for atmospheric stability throughout the HR- 
diagram. Note, however, that a sound analytical foundation for this numerical 'criterion' is 
currently lacking, and that important thermodynamic effects due to improved descriptions 
of Ti, as we present here, are crucial to ascertain its validity. Nevertheless, atmospheric 
regions where ionization and radiation cause I\ to locally decrease to below 4/3 are of great 
interest. These real gas effects also influence T 2 and T 3 in supergiants (e.g. Lobel et al. 1992). 
It is well-known that the local lowering of T 2 causes the onset of convective motions. Unsold 
(1948) mentions convective 'instability'-regions in the partial ionization zones of H and He. 
For dynamically stable atmospheres the local decrease of T 3 is linked with an increase or 
decrease of oscillation amplitudes over time, which determines the star's pulsational stability. 



7. Application to supergiant model atmospheres 
7.1. Model grid 

Based on the NLTE expression for Ti in Sect. 4.2 we investigate dynamic stability of 
atmospheric models we calculate for a range of effective temperatures and gravity acceler- 
ations for cool supergiants. The new model grid is computed to very high optical depths 
of log(rR OSS )~5 with ATLAS9 (R. Kurucz, private communication). A modified version of 
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ATLAS9 is used which, for certain models, includes 999 optical depth points. The number 
of points per decade has been improved for the computation of numerical derivatives in the 
treatment of convection. In deep model layers the convection decreases or completely van- 
ishes. More information on the ATLAS9 code and these plane-parallel hydrostatic models is 
available from Kurucz 1 web site. The temperature structure for models with spherical sym- 
metric geometry will change significantly, and more detailed solutions of the NLTE-problem 
(i.e. by solving detailed rate equations) will influence the behavior of IY In the present 
application we assume that the stellar radiation field in the optically thin part of the model 
atmospheres (tr oss <2/3) can be approximated by the stellar effective temperature for our 
computation of T 1 with height (see Sect. 8 for a discussion). 

The upper panel of Fig. 3 shows the behavior of F 1 (bold solid line) in the model with 
T cff =8000 K and log(g) = 1.0. In the partial NLTE-ionization region of hydrogen Ti decreases 
to ~0.8, around log(rR ss)=— 2. Towards smaller tr oss , r\ increases and assumes values of 
~4/3 (dotted horizontal line) in the outermost atmospheric layers. It results from the ionizing 
stellar radiation field which enhances the mean ionization fraction (solid line in the lower 
panel) of very optically thin layers. When only collisional ionization is considered (dash- 
dotted line), the local kinetic temperature becomes too low to appreciably partially ionize 
the fluid, and I\ exceeds unity. The Ti-minimum is therefore determined by the incident 
stellar radiation field which dilutes above R* (r Ross <2/3). The influence of lowering W(z) in 
the NLTE-ionization equation and in the radiation pressure is shown by the dashed lines for 
different values of z. The dilution of the ionizing radiation field diminishes the local value of 
x, and the point where x ~0.5 occurs in layers of increasingly smaller optical depth. Hence, 
the corresponding minimum in r\ (upper panel) displaces towards smaller TR OSS -values. 

At large optical depths (t Ross >2/3) the radiation temperature approaches the kinetic 
gas temperature. I\ decreases to below 4/3 due to thermal (LTE) ionization of hydrogen 
and helium, whereas a smaller decrease of I\ results from He + -ionization at the base of 
the model. At very large optical depths (r Ross > 10,000) we compute that the fluid becomes 
fully ionized, and I\ assumes constant values around 1.45. Deeper down at the base of the 
stellar envelope, the value of T± approaches asymptotically 4/3, because radiation pressure 
outweighs the gas pressure for high temperatures above 10 5 K. 

Figure 4 shows the behavior of Ti in model atmospheres for 4000 K< T e s <20,000 K, 
with small log(g)-values. The models are plane-parallel and assume a constant value of 
2 kms -1 for microturbulence with depth (Kurucz 1996). Although the effect of spherical 
geometry of extended atmospheres on the hydrostatic solution for the thermodynamic state 



1 http://kurucz/grids.html 
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is not negligible, we presently compare differences in T rad of at least 500 K, which are suf- 
ficient to outweigh this effect. We therefore infer global trends in r\ for a wide range of 
effective temperatures. In the model with T eff =4000 K (upper left panel) I\ assumes values 
around 5/3 in the outer atmospheric regions. In deeper layers Ti lowers to ~1.15, primar- 
ily due to thermal (LTE) ionization of hydrogen. For T e g=5000 K the stellar radiation 
field becomes sufficiently intense to partially photo-ionize hydrogen at very small optical 
depths, which strongly diminishes IY We set W=l/2 in our further calculations. When 
T eff is increased to 6000 K, 7000 K, and 8000 K (upper panel right), we find that the par- 
tial photo-ionization region displaces towards larger optical depths. I\ assumes very small 
local values ~0.8 for models with T rad =7000 K and 8000 K. This results from the thermal 
ionization zone of hydrogen, but which displaces towards smaller optical depths for models 
with higher T e fj. This combined effect whereby the partial thermal ionization zone displaces 
outwards, and the photo-ionization region displaces inwards by raising T e $, causes Y± to as- 
sume very small values (below unity) over a large geometrical fraction of supergiant models 
with 7000 K<T c{i <8000 K. 

In models with 9000 K< T c g <13,000 K (lower left panel) we find that the minima in 
Ti increase above i?^. This is because hydrogen becomes further ionized, or a;>0.5 towards 
higher T ra( j. For the model with T e g=ll,000 K, hydrogen becomes nearly fully ionized, and 
Ti exceeds 4/3 in the upper atmosphere. In these models the thermal He-ionization region 
occurs at increasingly smaller optical depths, and for T eff =13,000 K the region approaches 
r Ross =2/3 (~-R*). Around these effective (or radiation) temperatures, the partial photo- 
ionization zone of helium enhances in the upper atmospheric layers, and becomes noticeable 
by the decrease of Ti to below 4/3. The lower right-hand panel of Fig. 4 shows the model 
with T c ff=15,000 K, in which the thermal and photo-ionization regions of helium merge, 
which reduces r x to ~1.25 around R*. We find that atmospheric regions with Ti below 4/3 
are absent in models with T c r > 16,000 K. For the latter, the partial ionization of He + causes 
a minor decline in r 1; but which exceeds 4/3 throughout the entire model atmosphere. 



7.2. Model gravity dependence 

A comparison of Ti for models with extended atmospheres reveals that their dynamic 
stability, according Eq. (67), strongly depends on T e Q, or the ionizing stellar radiation 
field. However, differences in Ti with depth in these models also depend on the gravity 
acceleration. Figure 5 shows the contour-map of Ti in the (log(P g ), T c )-plane. We compute 
Ti with T rad =T cff , for models with 6000 K (upper panel) , 8000 K (middle panel), and 11,000 K 
(lower panel), by setting W—l/2. For T cff =6000 K, and P g between 0.1 and 1 dyn cm -2 , 



-22 - 



the ri-minima occur for 2000 K< T c <8000 K. This region marks the upper atmospheric 
layers where —6 < log(r Ross ) <— 4 in the model with log(g)=1.0 (see upper right-hand panel 
of Fig. 4). The solid diamonds in the upper panel of Fig. 5 show the (log(P g ), T c )-values 
of this model, whereas the filled triangles and boxes are plotted for models with log(g)=0.5 
and 0.0, respectively. The plots reveal that the revalues in layers of small optical depth 
and T e <4000 K, are not dependent of log(g), and assume comparable values determined 
by the local gas pressure. However, in layers with r Ross >2/3, where T e >T eff =6000 K, Ti 
assumes very different values ranging between 1.6 and 1.36, but which decrease towards 
smaller log(g)-values. The same trend is noticeable in the Ti contour-map for T e ff =8000 K. 
The solid diamonds are shown for log(g)=2.0, the triangles for 1.5, and the boxes for 1.0. 
Around tr oss =2/3 (where T e ~T e fj) Ti decreases to below unity towards smaller log(g)-values, 
as a result of enhanced partial photo-ionization of hydrogen. It also reveals that models 
computed for lower log (g) -values would assume even smaller revalues in the optically thin 
region of the atmosphere. However, for log(g)-values below 1.0 we could not converge models 
to a hydrostatically stable solution with T e g=8000 K. This results from the outwards directed 
radiation pressure which strongly increases the atmospheric density scale-height with smaller 
gravity acceleration, and which reduces <r x > to below 4/3. In the next section we show 
that the <Fi>-integral yields the atmospheric gravity values for which hydrostatic models 
of a given T eS become unstable. 

The lower panel of Fig. 5 shows the Tx contour-map for T eff =ll,000 K. The par- 
tial hydrogen ionization region, with small revalues, occurs at high kinetic pressures of 
log(P g )~1.5, because of the increased radiation temperature. The decrease of Tx in the par- 
tial ionization region of helium is also noticeable at these pressures around T c ~16,000 K. 
Three atmospheric models are plotted for log(g)=3.0 (diamonds), 2.5 (triangles), and 2.2 
(boxes). However, the latter has been extrapolated from the former two models. The de- 
tailed contour map reveals that the extrapolation does not lower Ti to appreciably smaller 
values (i.e. to below unity) in the atmosphere. It results from hydrogen becoming nearly 
fully ionized by the incident radiation field, whereby Ti assumes values around 4/3, also 
shown in the lower left panel of Fig. 4. 



7.3. Stability of model atmospheres 

We evaluate <r x > for our grid of model atmospheres. With dV = —dp/p 2 , Eq. (67) 
transforms to: 

<r ' >= ft-^h, ' (68) 

j pr p r 
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where /?r is the density of the outermost atmospheric layer, and p± is the density at the 
base of the stellar envelope. In layers of very high optical depth, beyond the partial H- and 
He-ionization zones, r\ becomes locally constant, and <r x > assumes an almost constant 
value. The models are dynamically stable against radial mechanic perturbations when <I\> 
exceeds 4/3. In general, we find that this condition is valid for models which we can converge 
to a stable solution. However, towards smaller log(g), <I\> can decrease to values very close 
to 4/3. 

Figure 6 shows the run of Ti with log(rR OSS ) (thin solid lines). The corresponding 
changes in <Ti> (bold solid lines) are shown by integrating over the atmosphere with Eq. 
(68), to the density at this optical depth. The atmospheric density structure of these cool 
supergiants is also shown for different log (g) -values of 0.5 (solid lines), 1.0 (short dash-dotted), 
1.5 (long dashed), and 2.0 (long dash-dotted). For these small gravity accelerations density 
inversions occur by the increase of continuum opacity in the hydrogen ionization region, 
which causes the outward directed radiative pressure gradient to exceed the gravitational 
pull. The condition of hydrostatic equilibrium therefore requires an outward increasing 
kinetic pressure (for a review see Maeder 1992). Pressure inversions occur in our models 
with 5250 K< T c g <8250 K, and 0.0< log(g) <1.5. The density-inversion layers occur at 
smaller optical depth with increasing T e fj, because the thermal hydrogen ionization region 
displaces outwards. 

Our calculations demonstrate that NLTE-ionization of hydrogen in the upper atmo- 
spheric layers of cool and extended atmospheres strongly diminishes <r x > to below 4/3. 
The value of <Ti> at a given depth point in Fig. 6 is a measure for the dynamic stability 
of the entire atmosphere above this point. A comparison of atmospheric conditions for dif- 
ferent T cff between 5500 K and 8500 K reveals that changes in the radiation temperature by 
1000 K strongly influence the behavior of <T X > with depth. The partial NLTE-ionization 
region of hydrogen displaces to higher depths by raising T eS from 5500 K to 6500 K (upper 
panels). The minimum of Ti moves deeper in the atmosphere, and consequently, <Ti> as- 
sumes smaller values around these depths. An integration of our model with T e ff=6500 K 
and log(g)=0.5 to higher depths yields a strong decrease of <r x > to below unity around 
log(r Ross )=1.5 (upper panel right). It partly results from the density-inversion which occurs 
around these depths. The stability integral Eq. (68) is determined by the local behavior 
of Ti and the atmospheric pressure- and density- structure. Towards larger optical depths 
<Ti> increases rapidly, because Ti assumes local values above 4/3 (dotted horizontal line), 
and the density increases steeply. Above r Ross ~1000, the stability integral assumes almost 
constant values of ~1.37 at the base of the stellar envelope. The integration reveals that 
this supergiant model, with very low gravity acceleration, is dynamically stable because the 
deeper evelope strongly contributes to the overall stability integral. However, we could not 
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converge a model with T er r=6500 K and a smaller log(g) of 0.3 to a stable solution. Models 
with very small gravity yield <r!>-values below 4/3 when integrating Eq. (68) to the base 
of the stellar envelope. 

Towards higher T er r of 7500 K (lower left panel) we compute that the minimum of r\ 
due to NLTE-ionization displaces further down the atmosphere, whereas the thermal H- and 
He-ionization regions occur at lower densities, higher in the atmosphere. However, in this 
model of higher T e g, but for the same log(g)=0.5, the gas density and pressure also diminish. 
Therefore, the inward integration of the local Ti causes a smaller decrease for <r\> at 
comparable optical depths, although I\ assumes locally very small minimum values of ~0.8. 
On the other hand, we find that by integrating this model of increased T ra d to very high 
depths at the base of the envelope yields a <rx>-value of 1.34, or only marginally above 
4/3. Stable models with log(g) below 0.5 could not be converged because <r x > falls to 
below 4/3 in the deepest layers. In general, we find that models with increased T eff (or 
^rad) become unstable towards higher gravity- values because <r\> drops to below 4/3. For 
T e fj=8000 K a stable solution is possible when only log(g)>1.0 (lower panel right), whereas 
for T e ff=5500 K we can converge stable models for log(g)>0.0. 



7.4. Atmospheric instability regions in cool supergiants 

Stable models of the same gravity yield towards higher T eS smaller <ri>-values at the 
base of the atmosphere, which is noticeable by comparing the four panels of Fig. 6. For a 
given T eS the stability integral increases with log(g), which demonstrates that the deeper 
regions of compacter models have a strong influence by stabilizing the overlying atmosphere. 
Towards very high log(g)-values the gravity pull effectively balances the radiation pressure 
gradient, and the density-inversion region vanishes. However, our NLTE-calculations of 
<Ti> reveal a remarkable property. <Ti> assumes very small values, around unity, down 
to the base of the atmosphere (up to r Ross ~100) for models with T eS around 6500 K — 7500 
K, practically independent of log(g). Models with higher or lower T eff (of the same gravity), 
yield larger <Fi>-values at these depths, or they have more stable atmospheres. 

To compare the run of Ti and <Ti> in models of different T e s and gravity, we plot the 
density scale in Fig. 7, rather than optical depths. The thin drawn lines show r 1; with the 
deep minimum due to NLTE-ionization of hydrogen, whereas the smaller minima at higher 
densities are due to the LTE-ionization zones. The curious loops in this scale result from 
the density-inversion for ln(p) between —24 and —20. The boldly drawn lines show the 
corresponding depression in <Ti> for models with log(g)=0.5 (solid lines) and log(g) = 1.0 
(broken lines). For the latter, <r x > assumes the smallest minimum of ~1.1 for ln(p)~— 24, 
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around T e g=6500 K. In contrast, the value of <r\> at the base of the model exceeds 4/3 
(m(p)> — 19), but it decreases steadily by further raising T eff in steps of 500 K. In other 
words, the base of these supergiant models becomes less stable towards T e g above 6500 K, 
whereas their extended atmospheric portions at lower densities tend to become more stable. 
The latter is also true for models with T e s below 6500 K, but instead the deeper envelope 
further stabilizes the entire model. 

Similar trends occur for models with log(g)=0.5. We have incorporated the stronger 
dilution of radiation in these more extended atmospheres by setting z—1.5, or the geometric 
dilution factor W=(3— V5)/6 (with Eq. 5) instead of 1/2, for the calculation of I\ above 
R*. The minimum in Ti therefore displaces towards smaller depths (see Sect. 7.1), be- 
cause partial ionization of hydrogen due to the more diluted radiation field occurs at lower 
densities. Consequently, <I\> assumes smaller values at lower densities. Our calculations 
reveal that for models with smaller gravity the minimum of <Ti> occurs for lower densi- 
ties, and vice versa. The dilution of radiation tends to further destabilize a more extended 
atmosphere above R+. In deeper layers where T ra d=T e , the atmosphere stabilizes because 
<Ti> increases, but assumes a smaller value above 4/3 (at the base of the envelope) than for 
models with higher gravity. Lower gravity models are less stable, and the regions which con- 
tribute to their destabilization occur at smaller densities in the atmosphere. These regions 
contribute strongest in supergiants with 6500 K< T eS <7500 K, as the result of partial non- 
LTE ionization of hydrogen with T rad ~T eff and the partial LTE-ionization of hydrogen and 
helium, combined with the complex mean density- and pressure-structure of such extended 
atmospheres. 



8. Discussion 

In a study of the evolution of luminous blue variable (LBV) stars Stothers & Chin 
(1994) evaluated <Ti> in the outer part of the stellar envelope and found that massive 
stars, evolving off the main sequence, develop dynamically unstable outer layers, for which 
eruptive mass-loss can result. Humphreys & Davidson (1994) pointed out that the low 
surface temperature of 12,000 K at which the first instability occurs in these calculations 
appears to be underestimated, because of the absence of very luminous late-B or A-type 
stars. We note that LTE calculations for our model with T e fj=12,000 K and log(g)=2.0 show 
indeed a decrease of r x to below 4/3 over a large fraction of the atmosphere with tr oss <2/3 
(see Fig. 5 of Lobel et al. 1992). Towards higher T e ff, hydrogen fully thermally ionizes, 
and hence Ti increases above 4/3. However, our present NLTE calculations of Ti show that 
stable supergiant models with T e s around 12,000 K have rather stabilizing outer layers (with 
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Ti slightly above 4/3; see Fig. 4) because the partial photo-ionization region where x~0.5, 
occurs deeper. 

We note, however, that we have set the radiation temperature equal to T eff , and W—l/2 
or (3 - y/E)/G in our calculations. More accurate computations of the detailed NLTE- 
ionization balance, in which T ra d varies with height over the atmosphere, and radiation 
pressure also dilutes with distance, are required to determine the precise boundary parame- 
ters where the models become dynamically unstable in the HR-diagram. Such calculations 
require 'case studies' of various individual supergiants. These stars can have very fast winds 
and high mass-loss rates, which strongly influences the atmospheric extension. Likewise, we 
expect that small deviations from an 'average' radiation temperature (but which is propor- 
tional to T eff ) will occur due to variations in the local opacity sources at different atmospheric 
levels. This T rad -dependency should be obtained from accurate observations of their spectral 
energy distributions. 

Our present calculations, which omit these further refinements, show however a clear 
trend in which <r\> assumes minimum values below 4/3 in stable models of supergiants 
with 6500 K < T eff < 7500 K. These atmospheres exist in the extension of the Cepheid 
instability strip, which has well-defined borders in the HR-diagram. In this area, the small- 
est values for <I\> occur down to the base of these atmospheres, practically independent 
of the gravity acceleration. It indicates the possible relation between pulsation-variability 
and dynamic destabilization mechanisms, caused by the decrease of <I\>. Soukup & Cox 
(1996) found with time-dependent calculations that the atmospheres of these massive stars 
become unstable to radial pulsations. They also mention that helium ionization dominates 
as the driving mechanism, and the contribution from the (thermal) hydrogen ionization zone 
becomes significant below 5000 K. Long-term spectroscopic and photometric observations 
of p Cas, with T cff =6500 K— 7250 K and log(L*/L©)=5.6— 5.9, demonstrate that non-radial 
pulsation modes are excited in these extended atmospheres (Lobel et al. 1994). In a theoret- 
ical study of oscillations in cool stars with log(L*/L )=5, Shibahashi & Osaki (1981) found 
that stable non-radial modes with Z=10 can be excited by the ionization zones. When the 
hydrogen and helium ionization zones are situated too shallowly they cannot excite radial 
pulsations, but they can drive non-radial modes, which are trapped near the stellar surface. 
For dynamically stable models, pulsation driving occurs in regions where T 3 approaches 
unity. Lobel et al. (1992) demonstrated that for conditions of cool supergiant atmospheres, 
the decrease of T 3 strongly couples with the decrease of Ti in the partial ionization zones. 
These considerations imply that the minima we compute for F 1 and <ri>, for low-gravity 
stars in the extension of the Cepheid strip, can be related with the excitation of stable non- 
radial pulsations by partial NLTE-ionization of hydrogen in the optically thin part of the 
atmosphere. 
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Furthermore, Shibahashi & Osaki (1981) mention convectively unstable zones located 
just below the photosphere due to hydrogen ionization in their model with T ef f=7080 K. In 
the model with T e fj=8000 K the ionization region emerges above the photosphere, which 
causes a strong difference between the modal properties of both models. In hydrodynamic 
simulations of compressible convection, which consider the effects of partial (LTE-) ionization 
of hydrogen, Rast (1992) mentions the importance of coupling of convection with pulsations 
in Cepheids. From an observational point of view, we note that photospheric absorption 
lines of yellow hypergiants display an unusually large macrobroadening, which cannot be 
attributed to rapid rotation (i.e. large wsim-values) for these evolved and very extended 
stars. For example, line profile modeling of high-dispersion observations in p Cas yield 
macro-broadening velocities of ~25 kms -1 (see Fig. 3 of Lobel et al. 1998). Fast large-scale 
velocities are possibly linked with ionization-induced formation of supersonic vertical flows 
(plumes), simulated for cool low-luminosity stars (e.g. Rast & Toomre 1993). 

Another remarkable aspect of yellow hypergiants (for a review see de Jager 1998) are 
'eruptions', which occur on time-scales much longer than the pulsation quasi-period (ca. 
half a century, say). In an outburst of 1945-46 p Cas (F8p) suddenly dimmed and displayed 
TiO-bands in its spectrum, characteristic of the photospheric temperatures of M-type stars. 
Within a couple of years (April '47) the star brightened up by nearly a magnitude, and 
a mid G-type spectrum was recovered around 1950. In 1985-86 the star showed a larger- 
than-average amplitude in the light curve, which could be associated with shell ejection 
events (Zsoldos & Percy 1991). Pulsational driving with the occurrence of strong convective 
motions may excite unstable modes with very fast growth rates, resulting in episodic mass 
ejection. If such mechanism can account for eruptions of yellow hypergiants, we expect that 
it is substantially different for the eruptions of LBVs, because strong convective motions 
induced by partial hydrogen ionization are not expected for these stars, and we compute 
<r 1 >-values above 4/3 over the entire atmosphere of hydrostatically stable models with 
T eff > 16,000 K. 

In a series of papers Nieuwenhuijzen & de Jager (1995), de Jager & Nieuwenhuijzen 
(1997), and Nieuwenhuijzen & de Jager (2000) presented strong theoretical and observational 
indications for the existence of a region in the upper HR-diagram where the atmosphere 
of yellow hypergiants become unstable. The cool boundary of the 'Yellow Evolutionary 
Void' occurs for bluewards evolving massive supergiants at T cff ~8300 K and log(L*/L©)>5.6. 
Evolutionary calculations show that these stars are expected to evolve along tracks of nearly 
constant luminosity, below the Humphreys-Davidson limit. They evolve bluewards because 
massive stars shed copious amounts of mass to the circumstellar /interstellar environment in 
the red supergiant phase. During redward evolution the very high mass-loss rates reduce the 
stability of the convective layer, and below a critical envelope mass it contracts into a thinner 
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radiative envelope which causes rapid blueward evolution. Only for a limited range of initial 
masses stars become red supergiants and evolve back far bluewards. A possible candidate for 
such a scenario is HR 8752, for which Israelian, Lobel, & Schmidt (1999) found an increase of 
the photospheric temperature by 3000 K— 4000 K, based on high-resolution spectra collected 
over the past 30 years. It is suggested that recurrent eruptions in yellow hypergiants occur 
when these stars approach the cool boundary of the Void, and 'bounce off' redwards. The 
bouncing against the Void may also explain why most of the cool luminous hypergiants 
cluster near its low-temperature boundary, while the identification of hypergiants of later 
spectral type, possibly like VY CMa (of M-type), is seldom. 

In their analysis of atmospheric acceleration mechanisms Nieuwenhuijzen & de Jager 
(1995) show that stable solutions cannot be computed for the atmospheres of blueward 
evolving yellow hypergiants at the low-temperature boundary of the Void. It results from 
a time-independent solution of the momentum equation, which considers the Newtonian 
gravity acceleration derived from the evolutionary mass, the gas, radiation and turbulent 
pressure gradients, and the momentum of the stellar wind. The latter requires the observed 
mass-loss rate. An 'effective acceleration' for the atmosphere is obtained iteratively, but 
which becomes negative at the cool border of the Void. The outwards directed net force 
causes an unstable atmosphere around T e fj=8300 K. We note that their momentum equation 
incorporates the important compressibility effects due to thermal ionization by the decrease 
of r\ in the adiabatic sound velocity, and of the isothermal and isochoric factors (x P and xt) 
in the gas pressure gradient. We suggest that the determination of instability boundaries for 
the Void be further improved with the NLTE-calculations of T 1 we present. It is important 
to note that in our computations the minimum of <I\> for stable models with lowest gravity 
acceleration occurs around T e ff=6500 K, which does not coincide with the low-temperature 
border of the Void. However, there is an overlap because <I\> approaches values very close 
to 4/3 in low-gravity models at this border. It suggests that the minima we compute for 
<I\> are rather linked with the driving of stable non-radial pulsation modes for the atmo- 
spheres of yellow supergiants in the extension of the Cepheid instability strip, whereas the 
Void can result from a particular combination of various acceleration mechanisms, combined 
with the decrease of the overall stability by the larger atmospheric compressibility due to 
partial NLTE- and LTE-ionization. The high mass-loss rates observed for these streaming 
atmospheres are also determined by the sonic point, which is situated inside the photospheres 
of cool massive supergiants. The dynamics of such extended atmospheres must therefore be 
linked with a prominent mechanism which drives these atmospheric pulsations, and which 
provides the momentum for their supersonic winds. 

Finally, we remark that more definitive conclusions about atmospheric variability mech- 
anisms should be obtained with fully time- dependent hydrodynamic simulations for a number 
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of prototypic hypergiants, like p Cas. Our analytical approximations are useful to evaluate 
important hydro dynamic quantities as the adiabatic sound speed = Ti Pt/p in conditions 
of NLTE, without having to numerically solve for large systems of non-linear balance equa- 
tions. With this analytical work we also provide a well-founded basis for more sophisticated 
numerical hydrodynamic calculations. Time-dependent codes which simulate the pulsation 
of extended atmospheres (i.e. described by Bessell, Scholz, & Wood 1996; Sasselov 1995) 
could adequately be updated for hydrodynamic NLTE effects, without the loss of their nu- 
merical efficiencies. With the present study we demonstrate that such advanced calculations 
cannot be based on the usual assumption of local thermodynamic equilibrium. In extended 
atmospheres, the thermodynamics of non-local equilibrium is required, and ultimately, these 
calculations also ought to consider non-equilibrium thermodynamic processes for realistic 
dynamic modeling. 



9. Conclusions 

1. We present new expressions for the first generalized adiabatic index I\ and the 
heat capacities, which are required for the study of dynamic stability of supergiant stellar 
atmospheres. Our equations consider important NLTE effects on the local ionization state 
by an incident and diluted radiation field, not in equilibrium with the local kinetic gas 
temperature. We demonstrate analytically that our more general expressions, which are 
also valid for multi-component gas mixtures, simplify to the classic formulae in which the 
radiation and kinetic temperature equilibrate. 

2. From a numerical application of our formalism to a grid of supergiant model atmo- 
spheres with solar abundance and 4000 K < T c r < 20,000 K, we find that the local values 
of Ti become very small, below 4/3, over a large fraction of the atmosphere for models 
with T eff between 7000 K and 8000 K. This results from the incident radiation field with 
a temperature of the order of the stellar effective temperature, which primarily causes par- 
tial photo-ionization of hydrogen. These regions displace deeper down the atmosphere with 
raising Teff, whereas the thermal partial ionization regions, at the base of the atmosphere, 
displace outwards. This combined effect causes the deep minimum in the local values of Ti to 
occur for these cool star atmospheres. Around these effective temperatures, partial NLTE- 
ionization causes very high atmospheric compressibilities, by which T 1 can even decrease to 
below unity. 

3. Numerical evaluations of Ledoux' stability integral <Ti> down into the stellar enve- 
lope demonstrate that <Ti> exceeds 4/3 in supergiant models for which we can compute a 
hydrostatically stable solution. Towards smaller gravity accelerations <Lx> decreases, and 
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models become unstable with <r\> < 4/3 at the base of the envelope. Stable models with 
smaller atmospheric gravity assume smaller <r!>-values, and Ledoux' stability integral ver- 
ifies that models of higher T eS destabilize at increasingly larger gravity-values due to the 
enhancement of radiation pressure. 

4. Most importantly, our calculations reveal that for stable supergiant models with 
6500 K < T c fj < 7500 K, <Ti> assumes minimum values, below 4/3, over a very large fraction 
of the atmosphere down to its base. Around T e g=6500 K this minimum occurs practically 
independently of the atmospheric gravity acceleration. These effective temperatures should 
be considered valid only within a case study for Kurucz atmospheric models. For a given 
T e ff, in atmospheres of increasingly smaller gravity, <I\> assumes small values in layers 
of increasingly lower density due to enhanced partial ionization with the dilution of the 
radiation field. This corresponds with dynamically destabilizing regions, extending over 
an increasingly larger geometric fraction in the upper atmospheres of more massive cool 
supergiants. 

I thank Prof. C. de Jager at SRON-Utrecht, for many discussions about the physics 
of supergiant atmospheres which have stimulated the development of this theoretical work 
over the past two years. Dr. R. Kurucz is gratefully acknowledged for calculating the new 
grid of atmospheric models and helpful discussions. The referee is thanked for several useful 
comments. This research is supported in part by an STScI grant GO-5409.02-93A to the 
Smithsonian Astrophysical Observatory. 
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Fig. 1. — The left-hand panels show the first adiabatic index I\ computed with partial 
NLTE-ionization for conditions of cool supergiant atmospheres. The mean ionization fraction 
x [panels right) varies with the radiation temperature and the kinetic gas temperature (T c ). 
Ti assumes minimum values for x~0.5., towards smaller T e . The upper panels are shown for 
a gas pressure of 10 dyn cm' 2 , and the lower panels for P g =0.5 dyn cm~ 2 . 
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Fig. 2. — Ti surface plots for P g =10 dyn cm 2 (upper panel) and P g =l dyn cm 2 (lower 
panel). The deep minima of F 1 result from partial NLTE-ionization of hydrogen with 
7000 K< T rad <9000 K. The minima between 14,000 K and 16,000 K result from the partial 
ionization of helium. The decrease of P g yields smaller minima for Ti with values below 0.7, 
which corresponds to very high gas compression ratios for these conditions in the atmospheres 
of cool supergiants. 
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Fig. 3.— The run of I\ in the model atmosphere with T efr =8000 K and log(#) = 1.0. The 
strong decrease of r x to below 4/3 (dotted horizontal line) results from partial NLTE- 
ionization for optical depths r Ross <2/3 (bold solid line). We assume T rad =T cff and 

W—l/2 for the computation of Ti above i?*. At larger depths in the photosphere, Ti mainly 
decreases due to thermal (LTE) ionization of H and He. The corresponding mean ioniza- 
tion fraction x is shown in the lower panel (bold solid line). At smaller optical depths the 
gas becomes fully ionized because T ra d exceeds T e . For unrealistic conditions of LTE, with 
T ra ,d=T e , these layers would assume a very small x, with larger revalues (dash-dotted lines). 
The influence of the dilution of radiation with distance z above R± is shown with dashed 
drawn lines. The mean ionization fraction assumes a;~0.5 at smaller optical depths when 
the ionizing radiation field dilutes more with increasing distance z. Consequently, the deep 
minimum of Ti in the upper panel, due to partial NLTE-ionization of hydrogen, occurs at 
smaller optical depths. 
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Fig. 4. — The behavior of T l in the outer envelope of supergiants with 4000 K< T cS <20,000 
K. Model T cS - and log(g)-values are labeled. The NLTE I\ assumes smallest values of ~0.8 
above i?* (tr oss <2/3) for models with T e g between 7000 K and 8000 K (upper panel right). 
We assume T ra d=T c g- and W=l/2 for the computation of r\ above i?*. Beneath R±, the 
LTE IV values decrease to below 4/3 (horizontal dotted lines) due to partial ionization of H 
and He. These thermal ionization zones displace outwards in models with higher T c g, while 
the partial hydrogen NLTE-ionization region in the outer atmosphere displaces inwards. 
This causes the deep ri-minimum for 7000 K< T cff <8000 K. Ti assumes increasingly larger 
values for models of higher T eS (lower panels). For models with T eff >16,000 K and log(g)>2.5 
helium is nearly fully ionized, and Ti assumes values above 4/3 over the entire atmosphere 
(see text). 
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Fig. 5. — Upper panel: Atmospheric models of supergiants with T e g=6000 K are plotted in 
the (log(Pg), T e )-plane, for three gravity accelerations of log(g) = 1.0 (diamonds) , 0.5 (trian- 
gles), and 0.0 (boxes). The solid lines show the contour map of Ti, computed with NLTE 
for T ra d=T c fj and W—l/2. The run of I\ in the outer atmospheric layers, with T c ~4000 
K— 5000 K, is very similar. The models differ more in the deeper layers, where different 
NLTE revalues are assumed. In layers with T e ~T efr (around R+), T 1 decreases towards 
smaller log(g)-values because the lines of equal Ti decrease parallel with the local (log(P g ), 
T c )-structure. In deeper layers (below R+), where T c >T cS , revalues evaluated with NLTE- 
ionization stay above 4/3. In these layers, LTE ionization calculations with T rad set equal to 
T c are required. 

Middle panel: Three models with T e ff=8000 K are plotted for log(g)=2.0 (diamonds), 1.5 
(triangles), and 1.0 (boxes). The NLTE Ti contour map is computed for T ra d=8000 K. In 
layers near R± (T e ~T e g), Ti assumes values below 4/3. The map reveals that in models with 
log(g) below 1.0 these layers will assume even smaller revalues, because their Ti-contours 
run parallel with the (log(P g ), T e )-lines. Stable hydrostatic solutions cannot be computed 
for these small gravity accelerations. 

Lower panel: Three models with T e gf=ll,000 K are plotted for log(g)=3.0 (diamonds) , 2.5 
(triangles), and 2.2 (boxes). In layers near R* (T c ~T cff ) the Ti contour map varies perpen- 
dicular to the atmospheric structure variations with gravity of the models. The layers above 
R± assume revalues around 4/3. Since hydrogen becomes nearly fully ionized, the upper 
atmospheres are more stable compared to the models with T e fj=6000 K and 8000 K. 
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Fig. 6. — The behavior of I\ (thin lines) and <I\> (bold solid lines) with optical depth in 
envelope models of cool supergiants for different gravity accelerations. <Ti> is integrated 
to the density at the corresponding optical depth. We assume T ra d=T c g and W—l/2 for the 
computation of r\ above R*. Higher gravity models are more stable because <r x > increases 
for a given T e g- at the base of the envelope with increasing log(g)=0.5 (solid lines), 1.0 (short 
dash-dotted), 1.5 (long dashed), and 2.0 (long dash-dotted). For a given \og(g), towards higher 
T eff , <ri> decreases at the base of the envelope to values close to 4/3 (dotted horizontal 
line). Models with log(#)<1.0 are unstable for T cff >8000 K. The model with T cff =8500 K 
and log(g) = 1.5 (short dashed lines in lower panel right) is also shown. The partial ionization 
zones of hydrogen and helium, where Ti locally decreases, are labeled. For T c g=6500 K and 
7500 K, <Ti> assumes minimum values down to the base of the atmosphere (log(rR OSS )~1.5) 
due to partial NLTE and LTE H- and He-ionization, combined with the density- and pressure- 
inversion regions of these models. Note the near gravity-independence for these minima in 
the models with T eff =6500 K (see text). 
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Fig. 7. — The run of F 1 (thin lines) and <r x > (bold lines) in the atmospheric density scale 
for log(g)=0.5 (solid lines) and 1.0 (broken lines). <Ti> assumes minimum values (below 
4/3) in the atmosphere of stable models with 6500 K< T cS <7500 K. We assume T rad =T eff 
for the computation of r\ above R*. Towards smaller gravity acceleration the destabilizing 
regions, due to partial NLTE-ionization of hydrogen, occur at lower density with the dilution 
of the ionizing radiation field (see text). 
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APPENDIX 

The first adiabatic index is defined by the adiabatic thermodynamic derivative: 
which can also be expressed by combining three thermodynamic quantities: 

ri = —x P , (2) 

c v 

where c v , cp t and xp denote the heat capacity at constant volume, at constant total pressure, 
and the isothermal factor, respectively. The total pressure, including the radiation pressure 
is given by the following equation of state: 

Pt = P g + ^rad = N k T c ( 1 + x ) p + ^ a W T r 4 ad , ( 3 ) 

where k is the Boltzmann constant, N the total number of particles per unit mass, and a is the 
radiation density constant. W is the geometric dilution function for the radiation pressure 
with distance d from the stellar surface R*. At the surface z—d/R^—l, hence W(z—\)—\/2 
with 

W(*) = 1/2(1- y/l-l/z*). (4) 

The ionization equation, modified for photo-ionization from the ground level for every ele- 
ment i, is (Sect. 2.1): 

i 

x = Q W ?L T rad exp (--j£A , (5) 

i-Xi p \ fcT rad y 

where Cj is a constant, and Xi the ionization fraction of element i, which is singly-ionized from 
the ground level, with ionization energy 7j, by the incident radiation field of temperature 
^rad- We denote the mean ionization fraction of the mixture with m elements: 



= ^2i/iXi, (6) 



x 

i=l 



where Vi = Ni/N is the element abundance having Ni particles per unit mass, and ^ i Vi — 1 
or N = J2i Hence, the total number of electrons per unit mass is N e = x N. 

When the radiation temperature T rad and the kinetic gas temperature T e are independent 
quantities, the specific heat capacity (per unit mass) at constant volume v (or density, 
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since v—l/p) is defined as the sum of the internal energy derivatives to both independent 
temperatures: 

' ' ■ (T) 



The internal energy function in Eq. (7) is given by: 



e = NkT e ^(l+x) + J2vi(. 



where x° = 1 — Xi denotes the neutral particle fraction. The two electronic energy terms 
of Eq. (8) include the partition functions which are defined by: Ui = 5^* e_ "^ an d 

Wi = J2T=o9r,i (fr") e _7 ^t , where g Tji is the statistical weight of the r th excitation level. The 
partition functions of the neutral fractions are indicated with the (0) superscript. 

In conditions of stellar atmospheres, the kinetic and ionic terms largely outweigh the 
electronic terms, and since we assume ionization from the ground level only, the electronic 
terms can be neglected, yielding: 

e = ?-NkT c (l +x) + N f>x t I»+ aWT ™ d . (9) 
1 i=i p 

The specific heat capacity at constant total pressure in Eq. (2) is defined as the sum of 
the derivatives of the enthalpy function h — e + P t / p to both independent temperatures: 

f dh\ ( dh \ 



x d T e J Pti Trad V d T rad / PuTc 
where m 4 4 

h=^-NkT e (l + x)+N Y j V i X i I i + 3 QH " r rad (n) 

We proceed by deriving the detailed expressions for the individual terms of Eq. (7) and 
Eq. (10). 

In Eq. (7) the kinetic temperature derivative at constant density yields: 

&X-H* + *<®) +K S«<&X- (12) 



and the radiation temperature derivative yields: 
\dT, 



— ^ = -NkT c (^-) +NYv i I i (4^ L ) +!2a9Nk(l + x), (13) 
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where a is the ratio of the radiation pressure and the kinetic gas pressure P ra d/-P, and 9 is 
the ratio of the local kinetic temperature and the temperature of the incident radiation field 

In Eq. (10) the kinetic temperature derivative at constant total pressure yields: 

dh\ 5 lr , . ^ f dx \ \ _ t ( dxj 

-Nk\(l + x)+T e 



dTj Pt 2 V \STjpJ j \9T,j ,, 

- iaNk{i+ ^{M) K (i4) 

and the radiation temperature derivative casts by using P ra( j — \aW T r ^ d into: 
dh \ 5 , r , „ / dx \ v-^ T / <9x,- 



C^rad / p t 2 y <9Prad /p. ^ \ C^rad / p t 

+4a ««( 1 + ,)(4-(^)J. ,15) 

Hence, both heat capacities are obtained from the detailed evaluation of two thermo- 
dynamic quantities in Eq. (12) and Eq. (13): 

9XA andfj^ , (16) 



dTj p ' \dT r3d/ r 



and with Eq. (6) 



(17) 

Four thermodynamic quantities are to be evaluated in Eq. (14) and Eq. (15): 

Oxi \ ( dxi \ f dp \ ( dp \ 

ydi) Pt ' {djtj Pt ' {dT e ) Pt ' and {df^ d ) Pt ' (18) 

for which also 



*9Te / p t ^ \ "9Pe / p t \ <9^~rad J Pt i V ^-^rad / p t 

i. The kinetic temperature derivative of the ionization equation (5) at constant density 
and radiation temperature yields after grouping terms and some rearrangement: 

9X< ^ *<!-*) I , P0) 



<91nT P / v M 2 \d\nT ( 



-51 - 



which casts with Eq. (17) (multiplying both sides with J2i v i )> after factorization, into: 

dx\ 1 x J2i - Xi) 



(21) 



We obtain by inserting Eq. (21) in Eq. (20) for every element v. 

dxA 1 Xi(l - Xi) Y^i v i x % 



0T c J p 2T C Y.i U i X i( l - X i) + Y,i V i X i' ^ 

Equation (12) reduces, after inserting Eq. (21) and Eq. (22) and grouping terms, to: 

+Nk ^^ i x i (l-x i ) 1 ^ r (23) 

i e 

= N k -(l + x) + -— / 2 " TaJ . (24) 

y2 K } 2 x + ^UiXiil-Xi) J 1 ; 

ii. Analogously, the radiation temperature derivative of the ionization equation (5) at 
constant density and kinetic temperature yields after grouping terms and some rearrange- 
ment: 

' =*<(l-*)((l + 7^V(^) ) (25) 



d\nT rad J p y\ kT rad J \d\nT iad/ p 

which casts with Eq. (17) after factorization into: 

dx \ 1 x Ei^i(l-^i) + 



dT rad J p T rad £\ u i x i( l ~Xi)+x 
We obtain by inserting Eq. (26) in Eq. (25) for every element v. 



(26) 



^ 1 x , {1 - x J( l + T ^]- ^- Xi) ^' \. (27) 



. d T iad J p T rad I V k T rad J Y.i "iZii 1 ~ x i) + J2i v i x 

Equation (13) reduces, after inserting Eq. (26) and Eq. (27) and grouping terms, to 

\dT Tad J p U2 V kT rad J x + ^UiXiil-Xi) 

+ JVfc^i/ i x i (l-xO^(l + ^)+12afl(l + x 



(28) 
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Hence, with Eq. (7) we obtain the total heat capacity at constant volume by summing Eq. 
(23) and Eq. (28), and grouping terms: 

= g + 12„«) (1 + x) + - «,) A g + .(i + 5 ^)) 



iVA; 



+ (^-£^(i-^J , + E,^d-^) ' (29) 

where we normalize the specific heat capacity to the gas constant to obtain dimensionless 
units. When denoting: 

X= ^^(1-2^), (30) 

i 

Y= ^^(l-a^A-, (31) 

i e 

and we define the functions Q and G for element i: 

1 

x + X ~~kT t 



Qi= i + ^fi + r^-). O 2 ) 



Ax-Y u , 
G * = ft + r^r » ( 33 ) 



Eq. (29) can formally be expressed: 

Q + 12 a fl) (1 + x) + "M 1 ~ x i) Gi ■ 



Nk 



(34) 



Equations (30) — (33) enable a similar compact form for the detailed expression of the heat 
capacity at constant total pressure, which we obtain below. 

Hi. The kinetic temperature derivative of the ionization equation (5) at constant total 
pressure and radiation temperature yields after grouping terms and some rearrangement: 

9X '' (35) 



d\nT e J Pt v ' \1 x(l +x) \d\nTjpJ ' 
which casts with Eq. (19) after factorization into: 

dx\ 3 x(l + x) Y^i^iO- ~ X i 



dT e J Pt 2T e x(l + x)+J2i^Ml-Xi) ' 



(36) 
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We obtain by inserting Eq. (36) in Eq. (35) for every element i: 

dxA 3 Xi(l - Xi) J2i v i x i (! + J2i u i x i) 



dT e J p 2T C J2i ^iXi{\ - Xi) + J2i u i x i (! + Ei v &i) ' 



(37) 



for which we note the formal resemblance with Eq. (22), apart from an extra factor (1 + 
Ei v i x i) m the numerator and the second term of the denominator. 

iv. The kinetic temperature derivative of the density at constant total pressure and 
radiation temperature results from differentiating the equation of state Eq. (3): 

9lnp ^ ^l + ^lMVl ) ■ (38) 



d In T C J Pt \ 1 + x \d\nT e/ Pt/ 

v. The radiation temperature derivative of the ionization equation (5) at constant total 
pressure and kinetic temperature yields after grouping terms and some rearrangement: 

dx- \ / /■ 1 / dx \ \ 

a^) P r xt{1 - x,) v + ^ + ^-w^)\9^)j • (39) 

which casts with Eq. (19) after factorization into: 

dx \ 1 x ( 1 + x ) ^iWiil-Xi) (l + 4a + 

. (40) 



dT rad J Pt T rad x(l + x) +J2iVi x i(l ~ x i) 

We obtain by inserting Eq. (40) in Eq. (29) for every element v. 

dx t \ = x) filial h Ei^(l-^)(l+4a+fcfc) 

d In T rad J Pt ' I k T rad J2i ViXiil - x^ + £\ v i x i i 1 + Ei 

(41) 

for which we note the formal resemblance with Eq. (27), apart from the extra terms with 4a 
in the numerator, and an extra factor (1 + £V i^Xj) in the second term of the denominator. 

vi. The radiation temperature derivative of the density at constant total pressure and 
kinetic temperature results from differentiating the equation of state Eq. (3): 

d\np \ / 1 / dx \ \ . 

1 '4a+- r — - . (42) 



91nT rad y Pt V. 1 + x \d\nT iad/ Pt/ 

Hence, with Eq. (10) we obtain the total heat capacity at constant total pressure by 
summing Eq. (14) and Eq. (15), and grouping terms: 
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5 m ( ( dx \ ( dx 



2 V V d T e J p \d 7r a( j j p 



k \\d T e y p t y<9 T rac j y p t 

where we normalize the specific heat capacity to the gas constant to obtain dimensionless 
units. 

When inserting Eqns. (36) — (38) and Eqns. (40) — (42) in Eq. (43), and defining with 
Eqns. (30) — (32) the function Hi for element i: 

we obtain, after grouping terms, the formal expression: 



iVA; 



Q + 4a (40(a + l) + l)^ (1 + s) + J] v iXi {l - Xi ) H t . 



(45) 



